/**
  ETFRoot project, Anyang Normal University && IMP-CAS
  \file cmd.C
  \brief a collection of commands for refined data analysis of ETFRoot
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2022-06-24
  \date 2021-06-24 last modified
  \attention
  changelog
  <table>
  <tr>  <th>Date         <th>Author      <th>Description                   </tr>
  <tr>  <td>2022-06-24   <td>Asia Sun    <td>file created                  </tr>
  </table>

  \copyright Copyright (c) 2021 Anyang Normal U. && IMP-CAS with LGPLv3 LICENSE
*/

const string dir = "/home/asia/Lab/ETFAna/build/";
const string file = "C15TA.root";
const string ff = dir+file;
const bool showCuts = false;

static const double PI = TMath::Pi();
static const string tstop0 = "tstop0"; // for cut C0
static const string TA = "TA"; // for cut B0
static const string pdcU1 = "pdcU1"; // for cutA0
static const double px[2] = {139.158, 0.40}; // A0 position and width x
static const double py[2] = {1.2184, 0.0423}; // A0 position and width y

void tmp(){
  TFile *f = new TFile(ff.c_str(), "UPDATE");
  TTree *ETF_TREE = (TTree *)f->Get("ETF_TREE");
  if(!ETF_TREE){ cerr << "nullptr ETF_TREE" << endl; exit(1); }
  //----------- define the A0, B0, C0 cuts ----------//
  const int n = 100; // number of points in a cut
  // CUT C0: boundary of the sensitive area of Tstop //
  TCutG *C0 = new TCutG("C0"); C0->SetTitle("C0");
  C0->SetLineColor(2); C0->SetLineWidth(2); C0->SetFillStyle(1000);
  C0->SetVarX("preTaTracker3D.GetHitPosX(tstop0)");
  C0->SetVarY("preTaTracker3D.GetHitPosY(tstop0)");
  C0->SetPoint(0,-27., 21.); C0->SetPoint(1,-27.,-25.); C0->SetPoint(2, 24.,-25.);
  C0->SetPoint(3, 24., 21.); C0->SetPoint(4,-27., 21.);
  // CUT B0: boundary of the target effective area //
  TCutG *B0 = new TCutG("B0"); B0->SetTitle("B0");
  B0->SetLineColor(2); B0->SetLineWidth(2); B0->SetFillStyle(1000);
  B0->SetVarX("preTaTracker3D.GetHitPosX(TA)");
  B0->SetVarY("preTaTracker3D.GetHitPosY(TA)");
  const double rTA = 22.;
  for(int i = 0; i <= n; i++){
    const double theta = 2.*PI*i/n;
    const double x = rTA*cos(theta) + 0., y = rTA*sin(theta) - 1.;
    B0->SetPoint(i, x, y);
  } // end for over i
  // CUT A0: choose the incident particle in dE0-tof1 spectrum
  const double xA0 = px[0], yA0 = py[0]; // ellipse position
  const double a = 3.*px[1], b = 3.*py[1]; // ellipse axis length: 3 sigma
  TCutG *A0 = new TCutG("A0"); A0->SetTitle("A0");
  A0->SetLineColor(2); A0->SetLineWidth(2); A0->SetFillStyle(1000);
  A0->SetVarX("pid.tof1()");
  A0->SetVarY("pid.dE0()+preTaTracker3D.GetHitPosY(pdcU1)*(0.005*pid.dE0()-0.0016)");
  for(int i = 0; i <= n; i++){
    const double theta = 2.*PI*i/n;
    const double x = xA0 + a*cos(theta), y = yA0 + b*sin(theta);
    A0->SetPoint(i, x, y);
  } // end for over i


  TCanvas *c = new TCanvas("c", "c", 800, 600); // backward check
  c->SetGrid();
  c->ToggleToolBar();
	c->ToggleEventStatus();
  if(showCuts){
    // confirm the cuts //
    ETF_TREE->Draw("preTaTracker3D.GetHitPosY(tstop0):preTaTracker3D.GetHitPosX(tstop0)>>(500, -60., 60., 500, -60., 60.)",
      "pid.dsca11()==0;tstop0X [mm];tstop0Y [mm];", "col");
    C0->Draw("same"); c->Print("C0.pdf");
    ETF_TREE->Draw("preTaTracker3D.GetHitPosY(TA):preTaTracker3D.GetHitPosX(TA)>>h(500, -60., 60., 500, -60., 60.)",
      "pid.dsca11()==0;TA_X [mm];TA_Y [mm];", "col");
    B0->Draw("same"); c->Print("B0.pdf");
    ETF_TREE->Draw("pid.dE0()+preTaTracker3D.GetHitPosY(pdcU1)*(0.005*pid.dE0()-0.0016):pid.tof1()>>h(500, 125., 150., 500, 0., 3.)",
      "pid.dsca11()==0&&C0;tof1 [mm];dE arb.;", "col");
    A0->Draw("same"); c->Print("A0.pdf");
  } // end if

  // return;

  // // --------------------- 1. MULTIPLICITY ------------------------ //
  // //// DC0-1-2 -- DCs downstream of the dipole magnet //////////////
  // ETF_TREE->Draw("Length$(dc0X0)>>h(13, -1.5, 11.5)");
  // ETF_TREE->Draw("Length$(dc0V1)>>h(13, -1.5, 11.5)");
  // //// PDCU0-1--D0-1 -- DCs around the target ////////////////
  // ETF_TREE->Draw("Length$(pdcU0X0)>>h(10, -1.5, 8.5)"); // pdcU0
  // ETF_TREE->Draw("Length$(pdcU1X0)>>h(10, -1.5, 8.5)"); // pdcU1
  // ETF_TREE->Draw("Length$(pdcD0X0)>>h(10, -1.5, 8.5)"); // pdcD0
  // ETF_TREE->Draw("Length$(pdcD1X0)>>h(10, -1.5, 8.5)"); // pdcD0
  // ETF_TREE->Draw("Length$(tofw)>>h(10, -1.5, 8.5)"); // tofw
  // ETF_TREE->Draw("Length$(preTaTracker3D)>>h(10, -1.5, 8.5)"); // track
  // // ------------- 2. FIRED WIRE/Strip DISTRIBUTION ---------------- //
  // ETF_TREE->Draw("dc0X0.GetWireId()>>h(83, -1.5, 81.5)"); // for DCs
  // ETF_TREE->Draw("dc0X1.GetWireId()>>h(83, -1.5, 81.5)");
  // ETF_TREE->Draw("pdcU0X0.GetWireId()>>h(10, -1.5, 8.5)"); // for PDCs
  // ETF_TREE->Draw("pdcU0X1.GetWireId()>>h(10, -1.5, 8.5)");
  // ETF_TREE->Draw("tofw.GetStripId()>>h(33, -1.5, 31.5)"); // for tofw

  // ------------- 3. Drift Time DISTRIBUTION ---------------- //
  // ETF_TREE->Draw("dc0X0.GetT()>>h(500, -100., 400.)"); // for DCs
  // ETF_TREE->Draw("dc0X1.GetT()>>h(500, -100., 400.)");
  // ETF_TREE->Draw("pdcU0X0.GetT()>>h(500, -100., 700.)"); // for PDCs
  // ETF_TREE->Draw("pdcU0X1.GetT()>>h(500, -100., 700.)");
  // // ------------- 4. Drift Distance DISTRIBUTION ---------------- //
  // ETF_TREE->Draw("dc0X0.GetR()>>h(500, -1.5, 6.5)");
  // ETF_TREE->Draw("pdcU0X0.GetR()>>h(500, -1.5, 11.5)");
  // // ------------- 5. Fired Channel DISTRIBUTION ---------------- //
  // ETF_TREE->Draw("tchpxi.GetChannelId()>>h(5001, -0.5, 4000.5)");
  // ETF_TREE->Draw("tchvme.GetChannelId()>>h(1501, 7499.5, 9000.5)");
  // // ------------- 6. TOT v.s. chid DISTRIBUTION ---------------- //
  // // ETF_TREE->Draw("tchvme.GetTOT():tchvme.GetChannelId()>>h(1501, 7499.5, 9000.5, 500, -100., 700.)",
  // // "tchvme.GetNofLeadingEdges()!=0&&tchvme.GetNofTrailingEdges()!=0");
  // ETF_TREE->Draw("tchvme.GetTOT():tchvme.GetChannelId()>>h(1501, 7499.5, 9000.5, 500, -100., 700.)", "tchvme.GetNofLeadingEdges()!=0&&tchvme.GetNofTrailingEdges()!=0", "col");
  // ETF_TREE->Draw("tchpxi.GetTOT():tchpxi.GetChannelId()>>h(5001, -0.5, 4000.5, 500, -100., 700.)", "tchpxi.GetNofLeadingEdges()!=0&&tchpxi.GetNofTrailingEdges()!=0", "col");
  // // Time To Trig
  // ETF_TREE->Draw("tchpxi.GetLeadingEdge():tchpxi.GetChannelId()>>h(5001, -0.5, 4000.5, 500, -100., 700.)", "tchpxi.GetNofLeadingEdges()!=0", "col");
  // ETF_TREE->Draw("tchpxi.GetTrailingEdge():tchpxi.GetChannelId()>>h(5001, -0.5, 4000.5, 500, -100., 700.)", "tchpxi.GetNofTrailingEdges()!=0", "col");
  // // ------------- 7. NofEdges v.s. chid DISTRIBUTION ---------------- //
  // ETF_TREE->Draw("tchpxi.GetNofLeadingEdges():tchpxi.GetChannelId()>>h(5001, -0.5, 4000.5, 8, -1.5, 6.5)", "", "col");
  // ETF_TREE->Draw("tchvme.GetNofLeadingEdges():tchvme.GetChannelId()>>h(1501, 7499.5, 9000.5, 8, -1.5, 6.5)", "", "col");
  // ETF_TREE->Draw("tchpxi.GetNofLeadingEdges()-tchpxi.GetNofTrailingEdges():tchpxi.GetChannelId()>>h(5001, -0.5, 4000.5, 9, -4.5, 4.5)", "", "col");
  // ETF_TREE->Draw("tchvme.GetNofLeadingEdges()-tchvme.GetNofTrailingEdges():tchvme.GetChannelId()>>h(1501, 7499.5, 9000.5, 9, -4.5, 4.5)", "", "col");
  // // ------------- 8. PXI Section Length DISTRIBUTION ---------------- //
  // ETF_TREE->Draw("PXISection.GetSectionLength()>>h(500, 0., 500.)");
  // // ------------- 9. Track Info ---------------- //
  // ETF_TREE->Draw("preTaTrackerX->GetK()>>h(500, -0.1, 0.1)");
  // // ------------- 10. Hit Position DISTRIBUTION ---------------- //
  // // pdcArrU
  // string d = "pdcU0"; ETF_TREE->Draw("preTaTracker3D.GetHitPosY(d):preTaTracker3D.GetHitPosX(d)>>h(500, -60., 60., 500, -60., 60.)", "", "col");
  // string d = "pdcU1"; ETF_TREE->Draw("preTaTracker3D.GetHitPosY(d):preTaTracker3D.GetHitPosX(d)>>h(500, -60., 60., 500, -60., 60.)", "", "col");
  // string d = "TA"; ETF_TREE->Draw("preTaTracker3D.GetHitPosY(d):preTaTracker3D.GetHitPosX(d)>>h(500, -60., 60., 500, -60., 60.)", "", "col");
  // string d = "tstop0"; ETF_TREE->Draw("preTaTracker3D.GetHitPosY(d):preTaTracker3D.GetHitPosX(d)>>h(500, -60., 60., 500, -60., 60.)", "", "col");
  // // pdcArrD
  // string d = "pdcD0"; ETF_TREE->Draw("postTaTracker3D.GetHitPosY(d):postTaTracker3D.GetHitPosX(d)>>h(500, -60., 60., 500, -60., 60.)", "", "col");
  // string d = "pdcD1"; ETF_TREE->Draw("postTaTracker3D.GetHitPosY(d):postTaTracker3D.GetHitPosX(d)>>h(500, -60., 60., 500, -60., 60.)", "", "col");
  // string d = "TA"; ETF_TREE->Draw("postTaTracker3D.GetHitPosY(d):postTaTracker3D.GetHitPosX(d)>>h(500, -60., 60., 500, -60., 60.)", "", "col");
  // string d = "tstop1"; ETF_TREE->Draw("postTaTracker3D.GetHitPosY(d):postTaTracker3D.GetHitPosX(d)>>h(500, -60., 60., 500, -60., 60.)", "", "col");
  // // postMagDCArr
  // string d = "dc0"; ETF_TREE->Draw("postMagTracker3D.GetHitPosY(d):postMagTracker3D.GetHitPosX(d)>>h(500, -400., 400., 500, -300., 300.)", "", "col");
  // string d = "dc1"; ETF_TREE->Draw("postMagTracker3D.GetHitPosY(d):postMagTracker3D.GetHitPosX(d)>>h(500, -400., 400., 500, -300., 300.)", "", "col");
  // string d = "dc2"; ETF_TREE->Draw("postMagTracker3D.GetHitPosY(d):postMagTracker3D.GetHitPosX(d)>>h(500, -400., 400., 500, -300., 300.)", "", "col");
  // string d = "tofw"; ETF_TREE->Draw("postMagTracker3D.GetHitPosY(d):postMagTracker3D.GetHitPosX(d)>>h(500, -400., 400., 500, -300., 300.)", "", "col");
  // // ------------- 11. PID ---------------- //
  // // incident particle identification //
  // // raw pid
  // ETF_TREE->Draw("pid.dE0():pid.tof1()>>h(500, 120., 150., 500., 0., 3.5)", "pid.aoz()>0.", "col");
  // // dE0 v.s. pdcU1 -- slant
  // string d = "pdcU1"; ETF_TREE->Draw("pid.dE0():preTaTracker3D.GetHitPosY(d)>>h(500, -50., 50., 500, 0., 3.)", "pid.dsca11()==0;pdcD0Y [mm];dE arb.;", "col");
  // // dE0 v.s. pdcU1 - corrected -- flat
  // string d = "pdcU1"; ETF_TREE->Draw("pid.dE0()+preTaTracker3D.GetHitPosY(d)*(0.005*pid.dE0()-0.0016):preTaTracker3D.GetHitPosY(d)>>h(500, -50., 50., 500, 0., 3.)", "pid.dsca11()==0;pdcD0Y [mm];dE arb.;", "col");
  // // dE0 v.s. tof1 - dE0 position dependence eliminated =====================================
  // string d = "pdcU1"; ETF_TREE->Draw("pid.dE0()+preTaTracker3D.GetHitPosY(d)*(0.005*pid.dE0()-0.0016):pid.tof1()>>h(500, 125., 150., 500, 0., 3.)", "pid.dsca11()==0;pdcD0Y [mm];dE arb.;", "col");
  // // confirm the cut - B0
  // string d = "TA"; ETF_TREE->Draw("preTaTracker3D.GetHitPosY(d):preTaTracker3D.GetHitPosX(d)>>h(500, -60., 60., 500, -60., 60.)", "", "col");
  // // confirm the cut - A0
  // string d = "pdcU1";
  // ETF_TREE->Draw("pid.dE0()+preTaTracker3D.GetHitPosY(d)*(0.005*pid.dE0()-0.0016):pid.tof1()>>h(500, 125., 150., 500, 0., 3.)", "pid.dsca11()==0&&C0;tof1 [mm];dE arb.;", "col");
  // A0->Draw();
  // // OFFTA and tofw slit - for B0 determination
  // ETF_TREE->Draw("pid.poz():pid.tof2()>>h(500, 60., 72., 500, 1300., 2200.)", "pid.dsca11()==0&&A0&&C0;tof2;poz", "col");
  // // &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& THE FINAL FRAGMENT PID SPECTRUM &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& //
  // ETF_TREE->Draw("pid.dE1():pid.aoz()>>h(500, 1.8, 3.7, 500, 0., 4.5)", "pid.aoz()>0&&pid.dsca11()==0&&A0&&B0&&C0;aoz;dE1 arb.", "col");
  // // ------------- 12. Tracking chi DISTRIBUTION ---------------- //
  // ETF_TREE->Draw("dc0X0->GetChi()>>h(500, -1., 1.)");
  // // STR
  // ETF_TREE->Draw("dc0X0->GetR()+dc0X0->GetChi():dc0X0->GetT()>>h(500, -100., 400., 500, -1., 6.)");
  // // chi-t
  // ETF_TREE->Draw("dc0X0->GetChi():dc0X0->GetT()>>h(500, -100., 400., 500, -1., 6.)");
  // // fit STR using exp data
  // h->ProfileX("pfx1"); pfx1->Fit("pol7", "R", "", 0., 350.);
  // // then collect the fitted parameters to the code
  //
  // // ------------- 13. Parallel Momentum Distribution ---------------- //
  // ETF_TREE->Draw("pid.dcospl()/pid.cospl()>>(500, -2.e-7, 5.e-6)");
  // // sigma(beta*gamma)/(beta*gamma)
  // ETF_TREE->Draw("2.00767e-3/pid.bg()>>(500, 0.002, 0.0036)", "A0");
  // // contribution of beta*gamma v.s. cospl to BGC (i.e. beta*gamma*cospl)
  // ETF_TREE->Draw("(2.00767e-3/pid.bg())/(pid.dcospl()/pid.cospl())>>(500, 10, 20000)");
  // // sigma(BGC)
  // ETF_TREE->Draw("pid.dBGC()>>(500, 1.868, 1.871)");
  // // sigma(BGC)/BGC: sqrt(sum2(DBG/bg, fdCospl/fCospl))
  // ETF_TREE->Draw("pid.dBGC()/pid.BGC()>>(500, 0.0022, 0.0034)");

  // ------------- 14. reaction products ------------------------ //
  ETF_TREE->Draw("pid.dE1():pid.aoz()>>h(500, 1.5, 3.7, 500, 0., 4.5)", "pid.aoz()>0&&pid.dsca11()==0&&A0&&B0&&C0;aoz;dE1 arb.", "col");

  const char *nucls[] = {"C15", "C15Ori", "C14", "C13", "C12", "B14", "B13", "B12",
    "B11", "B10", "Be11", "Be10", "Be9", "C15R", "C15U"};
  vector<TCutG *> cutls;
  for(const char *c : nucls) cutls.push_back((TCutG *)f->Get(c));
  for(TCutG *c : cutls) if(c) c->Draw("same");

  f->cd(); // return to the root folder

  return;

  // do the counting //
  for(TCutG *c : cutls) if(c){
    const double tmp = ETF_TREE->GetEntries((string("pid.dsca11()==0&&A0&&B0&&C0&&") + c->GetName()).c_str()); // &&(pid.dx2()>-9&&pid.dx2()<4)&&pid.dsca4()==1
		cout << c->GetName() << ": " << tmp << endl;
  } // end for over cuts
	cout << "N0: " << ETF_TREE->GetEntries("pid.dsca11()==0&&A0&&B0&&C0") << endl; // &&(pid.dx2()>-9&&pid.dx2()<4)&&pid.dsca4()==1

  // ------------- 14. reaction products -- END ------------------------ //








  // --------------- THE END ------------------- //
} // end member function cmd
